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Abstract 

Based on a PML for the advective wave equation, we propose two PML models for the 
linearized Euler equations. The derivation of the first model can be applied to other physical 
models. The second model was implemented. Numerical results are shown. 
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1 Introduction 



Since the work by Berenger on perfectly matched layer for the Maxwell equations Be r941 lBer96j 
in a computational box, many works have been devoted to a better understanding of their prin- 
ciple and behaviour see |MPV98| . [ZH96| . |CW94| . |LS00| |MC98j |BFJ 03| |BJ02] |A(;H02| to 
extensions to other geometries, see |ST04| |CM98j . or equations see |HN02j |A(;H99|[TTT03] . 
We consider here the linearized Euler equations which has been the subject of many works, see 
|Rahr)4| . |Hu01j . |TAC98j |Hes98j [TTu9fi] |Hag03| (and references therein). The key issue is a 



possible lack of long time stability [Hes98 , 



TAC98| . A first stable PML was proposed in jHuOlj 



for flows normal to the boundary. It has been recently extended to oblique flows in Hag03 



sec 
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also [FDFT02J. In these works, the PML is obtained via a change of coordinate in the complex 
plane applied to the direction normal to the boundary (say x). This amounts to replacing all the 
d x derivatives in the Euler system by an operator denoted by <9§ m ' which is still differential in 
the x direction but pseudo-differential in the other variables. Based on the Smith factorization 
WRLS2], we propose here another strategy for constructing stable PMLs. The main feature is 
that not all d x derivatives are modified. The Euler system is modified only in order to damp 
the modes that could reflect at the interface between the physical domain and the PML. In 
our approach, the vorticity modes which are convected by a transport operator are not damped 
since they leave naturally the computational domain at an outflow boundary. As a result, even 
after discretization, the vorticity in the PML region equals (up to the machine accuracy) the 
vorticity of the reference solution computed in a very large domain, see figure 01 in §0] 
Moreover, the technique introduced in this paper enables an appropriate "PML" treatment of 
the various scalar equations at the basis of the system of PDEs under consideration. We think 
therefore that the technique introduced in this paper should extend to other systems of PDEs 
as well (e.g. shallow water, anistropic elasticity, . . .). 

More precisely, in section [21 we introduce the Smith factorization of the Euler equation. This 
powerful tool is used in section 01 to propose two ways to design PML for the Euler equation 
(see DNR05 for the application of the Smith factorization in domain decomposition methods). 
Both PMLs are based on the use of a PML for the underlying advective wave equation. The 
derivation of the first model can be applied to other physical models. In section 01 numerical 
results are shown for the second model which is easier to implement. 

2 Analysis of the Euler system via Smith factorization 

We write the linearized Euler equations around a constant subsonic flow (u, v), a constant density 
p > and a constant speed of sound c > as: 

/ d t + ud x + vd y pc 2 d x pc 2 d y \ / P \ ( U \ 

l dx dt+ ud x +vdy \ \ U )= \ f U ) (1) 

\ \d y o d t + ud x + vd y J \V J \ f v J 

where {f p ,f u ,fv) T is a given right handside. 
2.1 Smith factorization 

We first recall the definition of the Smith factorization of a matrix with polynomial entries and 
apply it to systems of PDEs, see |Gan661 IGan98j (or WRL95 for the mathematical analysis of 
systems of PDEs) and references therein. 

Theorem 2.1 Let n be an integer and A an invertible n x n matrix with polynomial entries in 
one variable A : A = (aij(X))i<ij< n . 

Then, there exist three matrices with polynomial entries E, D and F with the following proper- 
ties: 

• det(E) and det(F) are real numbers. 

• D is a diagonal matrix. 
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A = EDF. 



Morevoer, D is uniquely defined up to a reordering and multiplication of each entry by a constant 
by a formula defined as follows. Let 1 < k < n, 

• S k is the set of all the submatrices of order k x k extracted from A. 

• Det k = {Det{B k )\B k £ S k } 

• LD k is a largest common divisor of the set of polynomials Det k . 
Then, 



Du(A) = TU^X}' l - k - n 



(2) 



(by convention, LDq = 1). 



Application to the Euler system We first take formally the Fourier transform of the system 
Q with respect to y and t (dual variables are k and lu resp.). We keep the partial derivatives 
in x since in the sequel we shall design a PML for a truncation of the domain in the x direction. 
We note 

iuj + ud x + ikv p<?d x ipcPk 



A 



Euler 



ik 
P 



iuj + ud x + ikv 

iu + ud x + ivk 



(3) 



We can perform a Smith factorization of A Euler by considering it as a matrix with polynomials 
in d x entries. We have 



where 



AEuler — EDF 

( 1 



D = 



and 



E 



(u(c z — u 



2))l/3 



1 

\ o o gt 

ipc 2 k 




(4) 
(5) 





u 



and 



F 



iuj + ud x + ivk E2 
X 






c 2 — u 2 



\ 



/ iui + ud x + ikv 

ikpc 2 ik 
o x iuj + ud x + ikv 



ikpc 2 ) 
l\ 



pu 
u 



V 



pu 



ioj + ikv 



iuj + ikv 



where 



Er> = u- 



-uc 2 + u 3 )d xx + (2u 2 — c 2 )(iuj + ikv)d x + u({iuj + ikv) 2 + k 2 c 2 ) 

c 2 (iuj + ikv) 



Q = iu + ud x + ikv 



(6) 
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and 

t = -J 2 + 2ikuvd x + 2iu(ud x + ikv) + (c 2 - v 2 )k 2 - (c 2 - u 2 )d xx 
The operators showing up in the diagonal matrix have a physical meaning: 

Q = d t + nS-j + vd y 

is a first order transport operator and 

C = d u + 2uvd xy + 2d t {ud x + vd y ) - (c 2 - v 2 )d yy - (c 2 - u 2 )d xx 
is the advective wave operator. 



(7) 



2.2 Modes via Smith factorization 

In the PML analysis of section we shall use the expression of solutions to the homogeneous 
Euler equation. In order to illustrate the previous section, we make use of the Smith factorization 
to compute them. We take the Fourier transform in t and y of (fT|) and seek non zero solutions 
to 

s / p(u,x,k) 

AEuier u(lo, X,k) 1=0 X G R, LJ G R, k G R 
\ v(u, x, k) 

Using Smith factorization @, we have 

p(uj, x, k) 
EDF | ii(uj,x,k) 
v(uj, x, k) 



xGR, wgR, /cgR 



Since det(-E) is a real number, E 1 is still a matrix with polynomials in d x entries so that we 
can apply it to the above equation and get: 



/ 1 
1 



\ o o gt 



This implies that 



p(oj, x, k) 
F | u(uj, x, k) 
v(u), x, k) 



p(u, x, k) 
F | u(uj, x, k) 
v(uj, x, k) 






Ei ai{uj,k)e 



xGR, uGR, fcGR 



Xi (ui,k)x 



x eR, wgR, fcGR 



(8) 



where gt{e x ^^ x ) = 0. Since QC is of third order in the x direction, we have three possible 
values for A^: 

X 1 = (9) 



u(iu + ikv) - c(iu + ikv)yjl - k ^ k ^P 
c 2 — u 2 



for \k\ \J c 2 — v 2 < \ui + kv\ 



(10) 



u(iu + ikv) — c-\Jk 2 (c 2 — v 2 ) — (u> + kv) 2 ) 



c 2 — u 2 



for \k\\/i 



v 2 > \u + kv\ 
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u(iu + ikv) + c(w + ikv) x /l - k ^%$ 

tor |fc| v c — f < |w + A;u| 



Aa = < 



c 2 — u 2 



(11) 



?/ 



{iuj + i/cu ) + cy / fc 2 (c 2 — w 2 ) — (a; + kv) 2 ) 
c 2 — u 2 



for \k\y c 2 — i> 2 > |cj + feu | 



Remark 1 Ai comes from the transport operator Q whereas A2,3 come from the advective wave 
operator C 

Since det(F) is a real number, i 7-1 is still a matrix with polynomials in d x entries so that we 
can apply it to equation (jHJ) and get: 

p(u,x,k) \ 3 / \ 

fi(w,ar,A:) J = ^ a«(w, k)F^ j | (12) 
v(ui,x,k) J 

We shall define, for % = 1,2,3 

Wi(w,Jb) = e-^C^F- 1 [ 1 (13) 




It is easy to check that indeed W% does not depend on x and that 



p(oo,x,k) \ 3 

fl(w,x,fc) ) = ^a i (w,A:)Ty j (^,A;)e Al( ^ )x (14) 
v(uj,x,k) J 



i=i 



In the classical mode analysis of a system of PDEs, vectors Wi(u, k) are obtained after the di- 
agonalization of a matrix. Using the Smith factorization simplifies the computation since these 
vectors are given by the explicit formula ([T5j) and not by an eigenvalue computation. 

Notation Let i and j be integers. For a matrix A, (A)ij denotes the entry of the i-th row and 
of the j-th column. For a vector V, (V)j denotes its i-th. component. 

Remark 2 In section W^l we shall use that 

0Oi3 = | (15) 

3 PMLs for the Euler System 

The Smith factorization of the Euler system Q and the computations of the previous section 
show that the modes correspond either to operator C or to operator Q. Among these two 
operators, the only operator which generates waves propagating in both positive x and negative 
x directions is the operator C. This suggests that designing a PML for the Euler equation can 
be reduced to the design of PML for the advective wave operator C 
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3.1 PML for the advective wave equation 

This question has been the subject of several works, |Hag03| , |H)HT02| . |HND2| . |BBBI)L04| 
|DJ03j |BBBDL03j and references therein. Following these works, we use the for operator C 
a PML defined by replacing the x derivatives by a "pml" x derivative. The definition is as 
follows. Let <t(cl>, x,k) > be the damping parameter of the PML and JF -1 is the inverse 
Fourier transform in the variables u and k, we define the pseudo-differential operator a{x) : 

/ ,,,, ^-1/ c(iu + ikv) . . 

a{x)(4>) = , V , 9 ^r—. 4>) 16 

wvr/ y c(iuj + ikv) + (c 2 -u 2 )a(u;,x,k) ' y ' 



We define the d x ml derivative by 

dl ml = a(x)[d x - ^^(d t + vdy)] + ^^(dt + vd y ) (17) 
c z — u A c z — u A 

We are now in position to define the PML-x equations of the advective wave operator £.: 

£ P mi = d tt + 2uvd y {dl ml ) + 2d t (ud p x ml + vd y ) - (c 2 - v 2 )d yy - (c 2 - u 2 m ml ) 2 (18) 
Let us notice that we have 

dr l ~d x = -y(x)[d x - ^^(d t + vdy)] (19) 
where the operator 7(3;) is a pseudo-differential operator in the t and y variables: 

tMM = ~t,\11 a{ ^ K ) H h (20) 

c{ilo + ikv) + (cr — u z )a{L0, x, k) 

A PML-y used for truncating the domain in the y direction would consist in replacing in the 
operator C the y derivatives by a "pml" derivative in the y direction defined as follows: 

dl ml = a(y)[d y - ^^(dt + ud x )} + ^z^ipt + ud x ) (21) 

The PML-y equations of the advective wave operator C for the truncation of the domain in the 
y direction read: 

£pml, y = dtt + 2uvdl ml {8 x ) + 2d t (ud x + vdl ml ) - (c 2 - v 2 ){dl ml ) 2 - (c 2 - u 2 )8 xx (22) 

In order to give a complete definition of a PML bordering a rectangular computational domain, 
we have three possibilities for the corner region. The first one consists in desiging a third PML 
model in the corner that is compatible with both PML- a; and PML-y as was done for the 
Maxwell system in Ber94 for instance. The second possibility is to use prismatoidal coordinates 
LS01 . The advantage is that it allows for arbitrary convex computational domains and not 
only rectangular ones. The third possibility consists simply in placing side by side PML-x and 
PML-y regions. We have implemented this last simple approach and obtain good results, see 
§0 Of course, the two first options deserve further investigations. 
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3.2 First PML model 

Based on (@J), a first possibility is to define a PML for the Euler system by substitution of C 
with C? ml in matrix D (see formula ©)• In matrices E and F and in the operator Q, the x 
derivatives are not modified. We modify only the advective wave operator. Let 



(23) 

\ o o gc v 

we define a PML for the Euler system 

= EDPmlp ( 24 ) 
with the following interface conditions between the Euler media and the PML 

C{{F{W E uler)h) = £pml((F(W pm l)) 3 ) 

Q{(F(W Euler )) 3 ) = G((F(W pml )) 3 ) (25) 
d x (G((F(W Euler )) 3 )) = (f x m \g{(F(W pml )) 3 )) 

where the subscript 3 denotes the third component of a vector. 

Study of the first PML media We now give the results of an analysis of the PML system 
similar to that of § 12.21 for the Euler system. We define 

\pml \ 

A 1 — Ai 



^,3 = C 1 + -/■ , ., -J (A2,3 - ^2 — + ivk)) + ^2 — - 



and for i = 1, 2, 3 

Wf ml \uj,k) = e ->T , ^)x F -i ( ] (26) 

\^ e \r l (^k) X j 

It is easy to check that indeed Wf ml1 does not depend on x and that for any solution W pml1 to 
the homogeneous first PML system there exist {j3i{uj,k))\<i<. 3 such that 

W (u,x,k) = ^2Pi(uj,k)W[ mll (io,k)e^ {uJ ' k)x (27) 
i=l 

If we consider the solution in the positive x half space, its boundedness as x tends to infinity 
implies that (3 3 = 0. 



3.3 PMLness of the first model 

A key property of a PML is that there is no reflection at the interface between the Euler media 
and the PML media. We will prove that it is the case for a truncation of the space with an 
infinite PML starting at x = with a x independant damping parameter a. We have to consider 
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the following coupled problem: 

Find (Wi,W r ) = ((pi,u h vi), (p r ,u r ,v r )) such that: 

A-EvlerWi = 0, t > 0, x < 0, y G R 

(28) 

ATuL W r = 0, t > 0, x > 0, y G R 
with the interface conditions (|25|) 

£((F(W Z )) 3 ) = £ pm K0W)) 3 ) 

a((F(iy / )) 3 ) = e((F(iy r )) 3 ) (29) 
9 2 (g((F(^)) 3 )) = dr\G{{F(w r )) z )) 

We take the Fourier transform in t and y of (|28jl and get: 

A E ulerWl = 0, x < 0, u, k G R 

(30) 

From section 12.21 we know that the general solution to the Euler system is : 

3 

Wi=J2 <*i(u, k)Wi(uj : k)e x '^ k > (31) 
i=l 

where W\ is defined in (|13fl . As for the solution in the PML media, we know from (|27|) and the 
boundedness of the solution as x tends to infinity that 

2 

ri r = Y,Pi(u,k)Wr ll (u,k)e x * ml ^ x (32) 
i=i 

We study the adequacy of the PML by considering a\ and ai to be given. This corresponds to 
ingoing waves from the Euler media and moving towards the interface between the Euler media 
and the PML media. The three other quantities (03, (/3^ ) i=1 i2 ) are determined by the interface 
conditions (|29|) . The media is perfectly matched if we have no reflection in the Euler media, i.e. 
if «3 = 0. We now prove that this is indeed the case. 
From pTj) and (|T3|) . we have: 

3 

(i ? (^)) 3 = E« l (^fc)e Al( "' fc)x (33) 

i=l 

From (|32j) and l)26|) . we have: 

2 

(F{W T )h = J2^^) eXrl{UJ ' k)X ( 34 ) 

i=l 

Let 

ciiio + ikv) 
a = 

c{iuj + ikv) + (c 2 — -u 2 )<r(u;, x, k) 
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Notice that 

Q = u(d x - Ai) 

t = _(c 2 - u 2 )(d x - \ 2 ){d x - A 3 ) 

C = -(c 2 - u 2 )a 2 (8 x - \l ml )(d x - Xf 11 ) 



and that 

so that the interface conditions yield: 



i ^ pml ^ pml 



Thus we have 



ai(Ai - A 2 )(A! - A 3 ) = ^(Af 1 ' 1 - ^i??* ~ 

a 2 (A 2 - AO + a 3 (A 3 - Ai) = (3 2 (X p 2 ml - Ai) (35) 
a 2 (A 2 - Ai)A 2 + a 3 (A 3 - Ai)A 3 = /? 2 (A^ m ' - Ai)A 2 

R 2 (Ai-A 2 )(A!-A 3 ) 

pi = a\a 



I \ pmi 1 \ pml \j\ pml \ pml \ 
A 2 — Ai 



02 = «2 

a 3 = 

The nullity of a 3 shows that there is no reflection at the interface between the Euler and an 
infinite PML media. 

In practice, the PML has a finite width. As a result, the coefficient /3 3 will not be zero in 
general. But, the corresponding mode W% e a x is exponentially decreasing as x is increasing. 
Its contribution to the solution on the interface can be made as small as necessary simply by 
increasing the width of the PML. 

Complexity of the first PML model A direct computation yields: 

/ 

ATJL = ^Euler + | (36) 

V Ci C 2 
where 

n (d x - cg ml )b[(u 2 - c 2 )(d x + dl ml ) + 2u(iu; + ivk)} , n C x 

L>\ = . — 9 , 7 . 1 ., — \ and u 2 = — 

ipc z k{iuj + ikv) pu 

The difference with the Euler system concerns only the last equation on the variable v, but : 

1. The formula is complex and involves third order derivatives on both the pressure p and 
the normal velocity u. 

2. The formula implies a division by ipc 2 k(iuj + ikv) which can be zero. 

As for the first point, one might argue that it is just a matter of implementation. The second 
point seems more serious. A possible cure could be to take: 



cf(uj, x, k) = &{x) (pc 2 k(iv + kv)Y 
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where a(x) > 0. From formulas for C\ and C2 and formula (|19j) - (|2()j) . we get then 



Ci = cr(x): 



(ipc 2 k(iui + ikv))(c 2 — v? 



-G[{u 2 -c 2 )(d x + di ml ) + 2 



U{IUJ 



ivk)] 



c(iuj + ikv) + (c 2 — u 2 )a(u>, x, k) ' 

As a result the symbols of C\ and C2 would not vanish. But it would be at the expense of the 
damping of the PML. Indeed, cr(uj, x, k) would be small for small values of k or of iu + ikv. The 
present first model raises difficulties. Nevertheless, it should deserve interest since it corresponds 
to a systematic way to design a PML for systems of PDEs. Moreover, since matrices E and 
F are not unique, it is quite possible that a more suitable Smith factorization when used in 
formula (|24|) would lead to a practicable PML. In the next section, we design another PML for 
the Euler system whose numerical results will be given in section |U 



3.4 Second PML model 

The rationale for this model is that the pressure p satisfies an advective wave equation which is 
the only equation that demands a PML. Indeed, let multiply © by the matrix 

I b -pc 2 d x 

El = 1 
\ 




(37) 



We get: 



El A 



Euler 



V 



p x 

ik 
P 





iuj + ud x + ikv 







ilo + ud x + ivk 



(38) 



We substitute £ with £ 



; pml 



and apply 



El 



-1 



V 



-pc 2 d x g- 1 
1 





and we are thus led to define: 



/ §-i(tp ml + c 2 (d xx -k 2 )) pc 2 d x ipc 2 k\ 



A 



pml2 
Euler 



p 

ik 
P 



Q 






h 



(39) 



A direct computation yields: 



Jprrt/2 
A Euler 



■^Euler + 



(t pml - b)g~ l 
00 
00 



In order to get rid of the operator Q \ we introduce a new variable V such that Q{V) 
that in the physical space the enlarged PML system we consider reads: 



(40) 



p so 



A' 



pml2 
Euler 



(V\ 

p 

u 

V v ) 



I 



Q 


-1 










£pml £ 


Q 


pc 2 d x 


pc 2 d y 





\d x 

p • L 


G 










-d 

p°y 





g 


) 



( V(t,x,y) \ 
p(t,x,y) 
u(t,x,y) 
\ v(t,x,y) J 



0, t>0,x>0,yGR 
(41) 
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with the following interface conditions between the Euler media and the PML 
V = 0, p and u are continuous, d x {pE u ier) = d x nnl {p pm i) 



Study of the second PML media We now proceed to an analysis of the PML system similar 
to that of § 12.21 for the Euler system. The Smith factorization of A p ^^ er reads 



Kir = EDF 



D 



(42) 



where 

10 
10 

o o g o 
V o o o gc,p ml j 

and E and F are matrices with polynomial in d x entries and their determinants are one. This 
will enable us to give the general form of the solutions to the homogeneous PML equations. 

Indeed, let us denote W = (P,p,u,v) T such a solution. From the Smith factorization, there 
exist (f3i(tjj, fc))i=o,...,3 such that 



F(W) = O 



( \ 


g A X 

\ / 



/ \ 

o 
o 



i=l 



V 



A? x 



e i 



where 



pml 


pml 
2,3 



pml 



Ai 
(1 + 



(c 2 - u 2 )a 
c{ioj + ikv) 



)(A 



2,3 



{iuj + ivk)) + 



;(itL> + ivk) 



By applying F 1 to the above equation, we see that there exist vectors W[ ml2 (uj, k), i = 0, . . . , 3 
such that 



ri = ^Tp t (u;,k)Wr l2 (u,k)e 



\? ml (ui,k)x 



(43) 



i=0 



If we consider the solution in the positive x half space, its boundedness as x tends to infinity 
implies that fa = 0. 



3.5 PMLness of the second model 

We proceed similarly to § 13.31 We will prove there is no reflection at the interface between the 
Euler media and the PML media for a truncation of the space with an infinite PML starting at 
x = with a constant damping parameter a. We consider the following coupled problem: 
Find (Wi,W r ) = {(p h u u v{), (V r ,p r , u r , v r )) such that: 

A Eu lerWi = 0, t > 0, x < 0, y G R 

A P Eul r Wr = 0, t > 0, x > 0, y G R (44) 
at x = 0, V r = 0, pi = p r , d x (pt) = #? mZ (p r ), ui = u r , t > 0, y G R 
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We take the Fourier transform in t and y of the above coupled system and get: 



A-EulerWl = 0, X < 0, LU, k € R 



%Zl r W r = 0, x > 0, lo, k e R (45) 

at a; = 0, P r = 0, ]5 Z = p n dzfj?/) = dx ml {p r ), in = u r , lo, k € R 
From section l2~2l we know that the general solution to the Euler system is : 

3 



W t = ^2a l (u J ,k)W l (u J ,k)e x ^ k > 
i=i 



where Wj is defined in (\YS\i. As for the solution in the PML media, we know from (|43|) and the 
boundedness of the solution as x tends to infinity that 

ri r = j2A(^k)Wr l2 (.u J ,k)e x r , (^ 

i=0 

We study the adequacy of the PML by considering a% and 02 to be given. This corresponds to 
ingoing waves from the Euler media and moving towards the interface between the Euler media 
and the PML media. The four other quantities (03, (/3i)i=o ... 2) are determined by the interface 
conditions. The media is perfectly matched if we have no reflection in the Euler media, i.e. if 
as = 0. We now prove that this is indeed the case. We focuse on the equation satisfied by the 
pressure. By applying matrix 



Q Q —pc 2 d x —ipc 2 k 
t° ^EuUri we nave that p r satisfies the equation of advective Helmholtz PML media: 

h pml {f> r ) = 

From (JHHJl) w e also have that 

hi>i) = 

The interface conditions on the pressure are pi = p T and d x (pi) = dx m \p r ). We know from 
works on PML for the convective Helmholtz that there is no reflection at the interface for the 
pressure. Therefore there exists j3 p {uj,k) such that 

^ = P p (u;,k)e x ^ k > (46) 

We prove now that as a consequence, is zero. Indeed, taking the first component of (fl3|) we 
have that 

3 

I* = J>i (W;)! (47) 

From (|13|) and 1)15(1 . we have 

(W i ) 1 = i^(e A ^) = i(A i -A 1 )e A ^ 
u u 
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So that we have (Wi)i = whereas (Wi)i 7^ for i = 2, 3. Thus we can infer from (|46|) and 
()47j) that 03 = 0. This shows that there is no reflection at the interface between the Euler and 
an infinite PML media. 

In practice, the PML has a finite width. As a result, the coefficient (3^ will not be zero in 
general. But, the corresponding mode W£ e s x is exponentially decreasing as x is increasing. 
Its contribution to the solution on the interface can be made as small as necessary simply by 
increasing the width of the PML. 

Remark 3 (PMLs for the 3D Euler system) Using obvious notations, we briefly mention 
that the design of a PML for the three-dimensional Euler system: 

( d t + ud x + vd y + wd z pc 2 d x pc 2 d y pc 2 d z \ 

jd x d t + ud x + vdy + wd z 

^d y d t + ud x + vd y + wd z 

y \d z d t + ud x + vdy + wd z J 

(48) 

is very similar to the two-dimensional case. Let G3 be the 3D first order transport operator and 
£.3 the 3D acoustic wave equation. The Smith form of the 3D compressible Euler equations is 



D = 



( 1 
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1 
















Gs 







V 








£3^3 


/ 



As in the 2D case, only the advective wave operator £3 needs a "pml" procedure. For instance, 
the second model should read: 
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<?3 
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pc 2 d y 


pc 2 d z 
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4 Numerical Results 

We have taken c = 300 and p = 1. The 2D linearized Euler equations are discretized on a 
uniform staggered grid using a Yee Scheme and a CFL equals to 0.3. The convective derivatives 
are discretized using an upwind scheme both in the Euler region and in the PMLs. The compu- 
tational domain is the square [0, 1.2] x [0, 1.2] and except in tabled PMLs have a width 5 = 0.9 
corresponding to n$ = 38 grid points. The damping parameter a depends on the coordinate 
normal to the interface (say x): o~(x) = a vm [ x 2 /<5 3 where a pm i is a positive constant. The initial 
solutions are zero. Let f(t,x,y) = (1 — 2i: 2 {f c t — l) 2 )e~' K ' 2 ^ ct ~ 1 ^ 2 8m{x, y) for t < T s and zero 
for t > T s with T s = 0.05, f c = A/T s and 5m is the Dirac mass located in the middle of the 
computational domain. The right handside was f(t,x,y) on all three equations of system (|T|) 
except for figurenwhere it was zero on the velocity components. The PML solution is compared 
with a reference solution that is computed on a much larger domain. In figure the pressure 
for both solutions are plotted as a function of time 4 points from the upperleft corner of the 
domain (left figure) and inside the PML (right figure). The velocity field is u = v = 1/3. In 
the Euler region, both curves are nearly identical. In the PML, we see the damping of the PML 
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50 100 150 200 250 300 50 100 150 200 250 300 

Figure 1: Pressure fields for the reference solution and the PML solution near the upperleft 
corner (left) and in the PML (right) for an oblique flow u = v = 0.33 vs. time steps 

solution. Of course, for the reference solution, this corresponds to an Euler region and there is 
no damping. In Figure 01 we have the same plot for the pressure and the vorticity as well, in 
this case u = 200 and v = 100. In agreement with the construction of the PML, we see that 
the vorticity is not damped at all in the PML region. The vorticity in the PML region equals 
that of the reference solution. Indeed, in the construction of the PML in (|38[) only the wave 
operator £ is "pml"-ized and the transport operator Q is not modified. In figure [5J we show 
the pressure at different times of the computation for an oblique velocity u = v = 270. For the 
same computation, pressure near the upperleft corner is shown on Figure0J FigureEJis a similar 
figure for a horizontal flow in a duct. 

In table d we study the influence of the parameters of the PML on the error between the 
reference solution and the pml solution. We see that with a generous PML (38 grid points), the 
error is indeed very small. With 18 grid points, the error is small and with 8 points results are 
not satisfactory even when using various parameters cr pm ;. It is worth noticing that due to the 
use of an upwind scheme for the transport operator, the error on the vorticity is equal to the 
machine accuracy for any layer parameters. 



Table 1: Relative errors in percentage for different PML parameters, (u = 200, v = 100) 



variable 


Relative error in percentage vs. (n$, cr pm i) 


(38, 40) 


(18, 40) 


(18, 80) 


(8, 40) 


(8, 80) 


(8, 160) 


P 


0.12 


1.0 


8.0 


10.0 


6.0 


1.0 


u 


0.25 


1.2 


0.4 


16.3 


0.4 


0.3 


V 


0.1 


0.1 


10.0 


20.8 


0.4 


10.0 


UJ 


2.10- 14 


3.10- 14 


2.10- 14 


2.10- 14 


4.10- 14 


2.10~ i4 



The long time stability of the PML was assessed by computing on time intervals five times 
longer than those used for generating the figures. No instability was observed for various flows. 
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Pressure 



Pressure 





Figure 2: Pressure field for an oblique velocity u = v = 0.9 at successive time steps 
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Pressures in the Euler region 
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400 
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800 
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Figure 3: Reference and "PML" solutions in the Euler and PML regions vs. time steps (u = 200, 
v = 100) 




Figure 4: Pressure field (left) and error on the pressure (right) near the upperleft corner for an 
oblique velocity u = v = 0.9 vs. time steps 
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Figure 5: Pressure field (left) and error on the pressure (right) near the upperleft corner for a 
horizontal flow (u = 0.33, v = 0) in a duct vs. time steps 

5 Conclusion 



The first PML model proposed in § 13. 21 is obtained by using the Smith factorization of the Euler 
equations and a PML for the advective wave equation. This method can be applied to other 
systems of partial differential equations (free-surface flow, anisotropic elasticity, . . .). The second 
PML model we have proposed for the Euler linearized equations are based on the PML for the 
advective wave equation. Thus, the PML for Euler inherits the properties from the latter. This 
second model was implemented and numerical results illustrate the efficiency of the approach. 
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